#### "Gender Inequality and Authoritarian Regimes" ####
# authors: "Lars Pelke"
# date: 2020-01-20

## use 01_data_wrangling before running this file 

## Set WD ##

# delete those countries with less than 3 country-year observations to dela with singular fit 
vdem <- vdem %>% # delete those countries with less than 3 country-year observations to dela with singular fit 
  dplyr::group_by(cown) %>%
  mutate(num_years = max(year) - min(year)) %>%
  filter(num_years >= 5) %>%
  ungroup(cown)

#### Figures 1 Main Paper ####

vdem$gwf_type <- as.factor(vdem$gwf_type)
vdem <- vdem %>%
  drop_na(gwf_type)

ggplot(data = vdem , aes(x = gwf_type, y =v2xpe_exlgender)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Political Exclusion by Gender")

ggsave("calculations/Output/Figure1.pdf", height = 14, width = 24, units= c("cm"))


#### Figure 2 Main Paper ####

ggplot(data = vdem , aes(x = gwf_type, y =v2xps_party)) +
  geom_violin(alpha=1, position = position_dodge(width = .75), size=1, color="black", fill = "grey90") +
  geom_boxplot(notch = TRUE, outlier.size = -1, color="black",lwd=1.2, alpha = 0.4, width = .2)+
  geom_jitter(shape = 16, size=1,  color="grey", alpha= 1) +
  theme_pubr() +
  labs(x = "Regime Types", 
       y = "Party Institutionalization")

ggsave("calculations/Output/Figure2.pdf", height = 14, width = 24, units= c("cm"))


#### Regime Types by Geddes et al and mulitparty competition ####

vdem <- vdem %>%
  group_by(cown) %>%
  fill(v2elmulpar_osp)
summary(vdem$v2elmulpar_osp) # compare L�hrmann et al. 2018 

vdem <- vdem %>%
  mutate(v2elmulpar_osp_bi = case_when(v2elmulpar_osp <= 1 ~ 0, 
                                       v2elmulpar_osp > 1 ~ 1, 
                                       is.na(v2elmulpar_osp) ~ 0))
table(vdem$v2elmulpar_osp_bi) # compare L�hrmann et al. 2018 

vdem <- vdem %>%
  mutate(gwf_type_mp = case_when(v2elmulpar_osp_bi == 1 & gwf_type == "party" ~ "party regime with ME", 
                                 v2elmulpar_osp_bi == 0 & gwf_type == "party" ~ "party regime", 
                                 v2elmulpar_osp_bi == 1 & gwf_type == "monarchy" ~ "monarchy with ME", 
                                 v2elmulpar_osp_bi == 0 & gwf_type == "monarchy" ~ "monarchy", 
                                 v2elmulpar_osp_bi == 1 & gwf_type == "personal" ~ "personal regime with ME", 
                                 v2elmulpar_osp_bi == 0 & gwf_type == "personal" ~ "personal regime", 
                                 v2elmulpar_osp_bi == 1 & gwf_type == "military" ~ "military regime with ME", 
                                 v2elmulpar_osp_bi == 0 & gwf_type == "military" ~ "military regime",
                                 v2elmulpar_osp_bi == 1 & gwf_type == "communist party" ~ "communist party with ME",
                                 v2elmulpar_osp_bi == 0 & gwf_type == "communist party" ~ "communist party")) %>%
  drop_na(gwf_type_mp)

table(vdem$gwf_type_mp)

vdem <- vdem %>%
  ungroup()

#### Prepare Data for analysis ####

vdem2 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, v2lgqugen, e_total_resources_income_pc, eprratio)

vdem3 <- vdem %>%
  drop_na(cown, year, gwf_type, e_wb_pop_ln, e_migdppcln, v2lgqugen, e_total_resources_income_pc, eprratio, v2xps_party)


#### Figures in the SUPPLEMENTARY APPEMDIX ####

vdem %>%
  group_by(year) %>%
  summarize(v2xpe_exlgender_mean = mean(v2xpe_exlgender), 
            v2xpe_exlgender_95 = quantile(v2xpe_exlgender, 0.95), 
            v2xpe_exlgender_05 = quantile(v2xpe_exlgender, 0.05)) %>%
  ggplot(aes(x = year)) +
  geom_linerange(aes(ymin = v2xpe_exlgender_05, ymax = v2xpe_exlgender_95)) +
  geom_line(aes(y = v2xpe_exlgender_mean)) +
  theme_pubr() +
  labs(x = "Year", 
       y = "Women's Political Exclusion")
ggsave("calculations/Output/FigureA1.pdf", height = 14, width = 24, units= c("cm"))

vdem %>%
  filter(year== 2010) %>%
  summarize(mean(v2xpe_exlgender))

#### Supplementary Appendix Descriptive Statistics ####

# Sample without Controls
cols1 <- c("v2xpe_exlgender", "gwf_type", "gwf_type_mp")

library(stargazer)
stargazer(
  title="Summary Statistics", 
  as.data.frame(vdem[, cols1]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"))

table_a1 <- vdem %>%
  group_by(gwf_type) %>%
  count() 

table_a1 <- vdem %>%
  group_by(gwf_type_mp) %>%
  count()
  

# Sample with Controls 
cols2 <- c("v2xpe_exlgender", "gwf_type", "gwf_type_mp", "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm",
           "e_wb_pop_ln_gm",  "eprratio_dm", "eprratio_gm", "e_total_resources_income_pc_dm", 
           "e_total_resources_income_pc_gm") 

library(stargazer)
stargazer(
  title="Summary Statistics", 
  as.data.frame(vdem2[, cols2]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"))


table_a2 <- vdem2 %>%
  group_by(gwf_type) %>%
  count() 

table_a2 <- vdem2 %>%
  group_by(gwf_type_mp) %>%
  count()


#Sample with Controls and Party Institutionalization 
cols3 <- c("v2xpe_exlgender", "gwf_type", "gwf_type_mp", "e_migdppcln_dm", "e_migdppcln_gm", "e_wb_pop_ln_dm",
           "e_wb_pop_ln_gm",  "eprratio_dm", "eprratio_gm", "e_total_resources_income_pc_dm", 
           "e_total_resources_income_pc_gm", "v2xps_party_dm", "v2xps_party_gm") 

# Model 5 and 6 #
library(stargazer)
stargazer(
  title="Summary Statistics", 
  as.data.frame(vdem3[, cols3]), type = "latex", 
  summary.stat = c("n", "mean", "median", "min", "max", "sd"))

table_a3 <- vdem3 %>%
  group_by(gwf_type) %>%
  count() 

table_a3 <- vdem3 %>%
  group_by(gwf_type_mp) %>%
  count()

rm(table_a1, table_a2, table_a3)

## Description Main Text ##
vdem %>%
  summarize(year = last(year))

vdem %>%
  count()


